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ABSTRACT 


The  Censored  Weibull  Program  (CENWEIB)  is  a UNIVAC  1108  computer 
program  for  fitting  censored  samples  to  a Weibull  distribution  or 
for  generating  random  data  from  a Weibull  distribution.  Maximum 
likelihood  estimating  equations  are  used  to  calculate  parameters 
for  the  two-parameter  Weibull  distribution.  Moment  estimating 
equations  are  used  to  calculate  parameters  for  the  three-parameter 
Weibull  distribution. 
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A COMPUTER  PROGRAM  FOR  FITTING  CENSORED  SAMPLES 
TO  A WEIBULL  DISTRIBUTION 
, (CENWEIB) 


1.  INTRODUCTION,  a.  The  Censored  Weibull  Program  (CENWEIB)*  is 
a computer  program  designed  for  use  on  the  UNIVAC  1108  computer. 

The  program  can  perform  either  of  two  specified  tasks.  First, 
given  the  scale  (a)  and  the  shape  ( B)  parameters,  it  can  generate 
data  from  a specified  Weibull  distribution.  Secondly,  it  can  fit  a 
Weibull  distribution  to  observed  data.  The  program  can  fit  data 
obtained  from  complete,  singly-censored,  or  progressively-censored 
samples.  Examples  of  data  are  the  failure  of  a life- test  compo- 
nent or  the  detection  of  a target  in  a time-to-detect  test.  Cen- 
sored data  are  obtained  when  units  are  removed  from  a test  before 
the  event  to  be  measured  (e.g.,  failure  or  detection)  has  oc- 
curred. There  are  two  types  of  censoring.  In  Type  I censoring, 
testing  is  terminated  after  a specified  time  interval.  In  Type  II 
censoring,  testing  is  terminated  after  the  specified  event  occurs 
a certain  number  of  times.  In  both  cases,  the  data  collected 
consists  of  the  termination  times  of  the  events  which  occurred, 
the  time  testing  is  terminated,  and  the  number  of  items  remaining 
when  testing  is  terminated.  When  all  remaining  tests  are  termi- 
nated at  one  specified  time,  the  test  sample  is  called  singly- 
censored.  Some  tests,  however,  are  conducted  so  that  at  the  in- 
itial censoring  only  a certain  number  of  items  are  terminated 
while  some  items  are  allowed  to  continue  until  the  designated 
event  occurs  or  until  they  are  terminated  at  subsequent  stages  of 
censoring.  Such  test  samples  are  called  progressively-censored. 
CENWEIB  can  be  used  for  either  Type  I or  Type  II  censoring, 
singly-censored  or  progressively-censored  samples,  as  well  as  on 
complete  (uncensored)  samples.  The  program  can  accommodate  up  to 
and  including  1,000  observations. 

b.  The  Weibull  distribution  can  approximate  a variety  of  dis- 
tributions. When  the  shape  parameter  (discussed  below)  of  the 
Weibull  distribution  is  1.0,  the  distribution  becomes  the  exponen- 
tial distribution.  When  the  shape  parameter  is  3.5,  the  distribu- 
tion closely  approximates  the  normal  distribution.  The  Weibull 
distribution  can  also  take  on  various  other  positively-  and 
negatively-skewed  forms.  CENWEIB  is  primarily  designed  for 
samples  which  are  believed  to  be  positively  skewed. 


*CENWEIB  was  initiated  by  Carl  B.  Bates  (CAA)  and  was  formu- 
lated and  programed  by  Keith  D.  Thorp,  a former  CAA  employee. 
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c.  The  basis  of  the  CENWEI8  Program  is  the  work  of  Cohen  (ref- 
erence 1)  and  Essenwanger  (reference  2).  Cohen's  scale  and  shape 
parameters  for  the  two-parameter  Weibull  distribution  are  calcu- 
lated by  using  maximum  likelihood  estimating  equations.  Simple 
iterative  techniques  may  be  used  to  derive  the  required  values. 

Once  Cohen's  parameters  are  found,  they  are  converted  to  the  more 
familiar  and  conventional  form  used  by  Essenwanger.  Essenwanger' s 

(moment  estimating  equations  are  used  to  calculate  the  parameters 

of  the  three-parameter  Weibull  distribution. 

d.  When  more  than  one  sample  of  data  is  generated  from  a 
Weibull  distribution  or  more  than  one  sample  of  data  is  fitted  to 
a Weibull  distribution,  the  CEMWEIB  Program  will  calculate  a Chi- 
square  statistic  which  can  be  used  to  test  the  statistical  equal- 
ity of  the  generated  or  fitted  distributions. 

2.  MATHEMATICAL  AND  STATISTICAL  DESCRIPTION,  a.  Weibull  Sample 
Generation.  CENWEIB  uses  the  Weibull  cumulative  distribution 
function  (cdf)  to  randomly  generate  Weibull  numbers  from  given  in- 
put parameters.  Uniform  random  numbers  are  generated  from  a uni- 
form random  number  generating  subroutine.  After  a uniform  random 
number  has  been  generated,  the  Weibull  cumulative  distribution 
function 

F(x)  = /q  b t^"^exp(-t/a)  dt/a®,  [1] 

is  integrated  from  0 to  x,  where  x is  the  point  on  the  distribu- 
tion that  defines  an  area  under  the  curve  equal  to  the  value  of 
the  uniform  random  number  which  was  generated.  N of  these  Weibull 
values  are  generated.  Then  using  these  N values  CENWEIB  calcu- 
lates estimates  (a  and  Bl^of  the  parameters  (a  and  B)  of  equation 
[1].  Substituting  a and  S for  a and  B,  gives  an  estimate  of  the 
probability  density  function  (pdf) 

i 

f(x)  = (e/a^)x^-lexp[-(x/a)®],  x > 0, a > 0,  B > 0.  [2] 

t 

, Random  variation  may  actually  cause  the  true  Weibull  distribution 

parameters  produced  by  the  generated  data  to  be  different  from 
those  originally  inputted. 

b.  Cohen's  Maximum  Likelihood  Estimation 

(1)  Parameter  Estimation.  The  method  used  to  estimate 
Weibull  parameters  in  the  development  of  the  two-parameter  Weibull 
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distribution  in  CENWEIB  is  Cohen's  maximum  likelihood  (ML)  estima- 
tion technique.  In  the  development  of  this  technique,  Cohen  uses 
the  following  forms  of  the  Wei  bull  pdf  and  cdf: 


f(x)  = (Y/e)x^"^exp(-  Ve);  x > o,y>  o,e>  0, 
F(x)  = 1 - exp(-x^/0), 


where  6 is  Cohen's  scale  parameter  and 
Y is  Cohen's  shape  parameter. 

Cohen  develops  his  parameter  estimating  equations  from  the  maximum 
likelihood  function.  For  complete  samples  the  likelihood  function 
is: 


L(x,,...,x  ;y,6)  = n (Y/0)Xi^exp(-xy/9) , 

i=l 


where  Xj^,  x2,....,Xn  are  all  uncensored  observations.  Taking  the 
partial  derivatives  of  the  natural  logarithm  of  equation  [5]  with 
respect  to  Y and  6 and  setting  the  results  equal  to  zero,  we  get: 

L n n in 

-- 2;x^ln  x^  = 0,  [6] 

^ ^ i=l  ” i=l 

In  L n 1 " Y 

-r— = --  + -2  = 0.  [7] 


8 6 


Eliminating  6 between  equations  [6]  and  [7]  we  obtain: 


n Y, 

1 xjln  Xj 


.-^/i'"  ^i  1 . 1 

1=1 - — Xln  X,. 

n 1 


I Xi 
1=1 


" “l 


Using  standard  iterative  procedures,  this  can  be  solved  for  the  ML 
estimate  y.  After  solving  for  y,  0 can  be  determined  using  equa- 
tions [6]  and  [7]  so  that 


i 

[ 
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n Y 


[9] 


The  denotes  an  estimator.  Maximum  likelihood  estimating  equa- 
tions are  analogous  for  Type  I and  Type  II  censoring.  For  singly- 
censored  samples,  the  ML  function  is: 


L = 


x|"^exp^--^jj[l  - F(xj)]^’", 


[10] 


where  N is  the  total  number  of  observations  and 

n is  the  total  number  of  uncensored  observations. 


F(xj)  is  the  Wei  bull  cdf  at  the  termination  point,  x-j-.  Then, 


31n  L _ n 

“3y  Y 


n 

Sin  x^ 
i=l 


^ Y 

S X^ln  x^  = 0, 


[11] 


91n  L 
30 


[12] 


where  X*  denotes  a summation  over  the  whole  sample  with  the  cen- 
sored observations  being  assigned  the  value  Xj. 


From  these  equations,  we  get 

SVln  Xi  1 1 

— "17 — ~ — = — Sin  x^ . 

I**;  V " ,7,  ' 

Again,  Y can  be  calculated  using  iterative  techniques,  and 


[13] 


[14] 


For  Type  I progressively-censored  samples, 

n k 

L = C n[l  - P(T^)]  i [15] 

7=1  i=l 
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where  C is  a constant, 

k is  the  number  of  times  censoring  occurred, 


(i  = are  the  times  of  censoring,  and 

r^  is.  the  number  of  observations  randomly  terminated  at  the 


i^”  stage  of  censoring.  Then, 


91n 


=7^.  II" 

^ i=l 


= 0, 


[16] 


9ln  L 
96 


n 1 XT’'*'*'  Y 


e 9 


[17] 


where  ^**  denotes  summation  over  all  the  observations,  with  an  ob- 
servation censored  at  time,  T.j,  assigned  the  value  T.j. 

We  can  then  derive: 

X, 

i**x; 

After  determining  y, 

0 = I**x|/n.  [19] 


1 1 y , 

1- 


[18] 


The  two  estimating  equations  for  the  Type  II  progressively- 
censored  samples  are  analogous  to  equations  [18]  and  [19]  although 
intermediate  steps  in  their  derivation  differ.  As  can  be  seen 
from  the  above,  the  likelihood  equations  for  each  of  the  three 
censoring  cases  (complete,  singly-,  and  multiply-censored)  are 
different.  However,  the  parameter  estimating  equations  derived 
from  maximum  likelihood  equations  are  basically  the  same.  This 
enables  us  to  use  one  equation,  each  with  the  appropriate  summa- 
tion, to  estimate  the  scale  and  shape  parameters. 

(2)  Variance-Covariance  Matrix.  Cohen's  variance-covariance 
matrix  is  approximated  by  using  the  estimated  values  of  the  pa- 
rameters to  construct  the  information  matrix.  The  inverse  of  the 
information  matrix  is  the  variance-covariance  matrix.  The  ap- 
proximate variance-covariance  matrix  is  thus: 
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9^1  n L 

9^1  n L 

-1 

’ V(Y) 

A A -■ 

cov(Y,e) 

92ln  L 

/V  A. 

Y,e 

9y30 

92in  L 

<CD 

< 

= 

Cov(Y,9) 

V{0) 

9e9Y 

A A 

Y,9 

902 

A A 

Y.e 

Calculation  of  the  second  parti als  of  the  complete  sample  ML  func- 
tion for  the  information  matrix  gives: 


9^1  n L 


n n /V 

9.6 

i=l 


9^1  n L 


1 " " 

= I x|ln  x^., 


9y99  Y.e  0 1=1 


9^1  n L 


9e‘ 


JL  i_  V Y 

C S = ■ -2  + -3  - >^i  • 

0 0 i=l 


[21] 

[22] 

[23] 


The  second  partial s of  the  ML  functions  for  the  singly-  and 
progressively-censored  samples  can  be  taken  in  a similar  manner  to 
obtain  results  analogous  to  equations  [21],  [22]  and  [23]. 

c.  Essenwanger' s Moment  Estimation 


(1)  Parameter  Estimation.  In  Essenwanger' s more  conven- 
tional notation,  the  Meibull  pdf  and  cdf  are: 


f(x)  = (6/a^)x®-lexp[-(x/a)^];  x > 0,a>  0,  6 > 0 [24] 

F(x)  = 1 - exp(-x/a)®  [25] 


where  a is  Essenwanger' s scale  parameter  and 
3 is  Essenwanger' s shape  parameter. 


CENWEIB  calculates  Wei  bull  parameters  based  on  Cohen's  form  of  the 
Weibull  distribution  and  converts  to  Essenwanger's  form  using  the 
relationships  3=  y and  a = 9^'^.  In  other  words,  Cohen's  and 
Essenwanger's  shape  parameters  are  the  same,  but  the  scale 
parameters  are  different,  although  related  to  each  other.  CENWEIB 
outputs  Essenwanger's  parameters. 
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(2)  Variance-Covariance  Matrix.  Essenwanger' s variance- 
covariance  matrix  is: 


(-a^in  -i)cii  + 


(4--  1) 


(-a^ln  -)cii  + — ^Ci2  [(-a^ln  — yci2][-a^ln  y] 


+ [(-cx^ln  l)Ci2  +(-a!i^jc22]Ca!Lj^] 

[26]^ 


(4  - 1) 


where  the  c's  are  the  elements  of  Cohen's  variance-covariance 
matrix: 


Cov(9.e) 


Cov(Y,0) 


(3)  Distribution  Mean  and  Variance.  The  mean  and  variance 
of  the  distribution,  in  Essenwanger's  notation,  are: 


y = aPd 


where 


= a^ini  +4^  - [HI  +4)]^), 


rix)  = f t^“^exp(-t)dt 


for  a given  value  of  the  argument  x. 
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(4)  Three-parameter  Wei  bull.  A segment  of  CENWEIB  calculates 
parameters  for  a three-parameter  Weibull  distribution. 
Essenwanger's  form  of  the  three-parameter  Weibull  distribution 
function  is 


F(x)  = 1 


exp[-(-^)  J, 


where  6 is  the  location  parameter.  The  location  parameter  identi- 
fies the  starting  point  of  the  function  on  the  abscissa.  The  two- 
parameter  Weibull  form  will  force  the  starting  point  of  the  func- 
tion to  be  at  zero  on  the  abscissa.  Since  the  three-parameter 
form  does  not  force  this,  it  allows  the  fitting  of  a curve  to  a 
set  of  data  which  may  more  closely  approximate  the  data  values 
than  the  curve  produced  from  the  two-parameter  form.  Although 
CENWEIB  plots  only  the  two-parameter  Weibull  distribution,  the 
parameters  of  the  three-parameter  Weibull  distribution  are  also  in 
the  output  as  an  aid  to  the  user.  Estimates  of  the  three  parame- 
ters are  found  using  Essenwanger's  moment  technique.  The  shape 
parameter  moment  estimating  equation  is 

c - Sab  + 2a^ 


(b  - a2)3/2  ’ 


where  a = Fd  + -g) 
D =r(l 
c = Fd  + |) 


As  can  be  seen  in  the  above  equation,  6,  the  shape  parameter,  is 
the  only  unknown.  It  can  be  found  by  employing  iterative  tech- 
niques. Once  3 has  been  found,  the  scale  parameter  is  given  as 

a =Va^/(b  - a^) 


and  the  location  parameter  is  given  as 


6 = y - a a 


where  y is  the  population  mean. 
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d.  Hypothesis  Testing.  When  more  than  one  sample  of  data  is 
generated  or  more  than  one  sample  of  data  is  fitted  to  a Wei  bull 
distribution,  CENWEIB  calculates  a Chi-square  statistic  with  two 
degrees  of  freedom  in  order  to  test  the  statistical  equality  of 
the  distributions.  To  compare  the  two  distributions,  the  null 
hypothesis 


- 

Hq: 

“l 

“2 

ca 

1 

32 

is  tested  against  the  alternative  hypothesis 


[32] 


- 

* 

“1 

°‘2 

^1 

32 

_ 

- 

The  hypothesis  is  tested  by  calculating 


[33] 


Q = [a^  - 02,  3^  - $2] 
where  the  variance-covariance  matrix  is 


a2(S)  ala,B)  ' 

-1 

A A 

oji  - 02 

a(a,3) 

3i  - 32 

[34] 


- " 

** 

a2(a)  0(0,6) 

V(aj^)  + V(a2) 

Cov(a^,6j)  + Cov(a2,32) 

, A A , 0 

0(0,6)  0^(6) 

Cov(aj^,6j) 

V(6i)  + V(62) 

+ Cov(a2,62) 

.[35] 


The  Q-statistic  is  approximately  distributed  as  a Chi-square  vari- 
ate with  two  degrees  of  freedom,  i.e.,  Q ~ (see  references  3 

and  4 for  more  detail).  An  inspection  of  the  Q-statistic  shows 
that  close  agreement  between  the  two  distributions  yields  a small 
statistic,  while  a large  difference  between  the  two  yields  a large 
statistic.  Therefore,  to  test  the  null  hypothesis,  compare  Q with 
X‘^(l-a,  2).  If  Q >x^(l-a,  2),  reject  the  null  hypothesis  at  the 
a-level  of  significance;  otherwise,  do  not  reject  the  null  hy- 
pothesis. By  rejecting  the  null  hypothesis,  we  are  saying  that 
the  two  distributions  are  not  equal. 
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3.  COMPUTATIONAL  PROCEDURE,  a.  Appendix  D contains  the  flow- 
chart for  the  program  and  a complete  listing  of  the  computer  pro- 
gram. First,  the  program  checks  to  see  if  data  are  to  be  fitted 
or  if  data  are  to  be  generated.  If  data  are  to  be  generated,  the 
program  generates  the  required  data  with  specified  scale  and  shape 
parameters,  orders  the  data,  and  censors  it.  If  data  are  to  be 
fitted,  the  program  first  orders  the  data,  then  calls  the  WEIBUL 
subroutine.  Using  the  data,  this  subroutine  calculates  parameters 
for  the  three-parameter  Wei  bull  distribution  using  moment  estimat- 
ing equations.  An  iterative  procedure  is  used  to  calculate  the 
parameters  if  there  are  censored  data.  Then,  Cohen's  parameters 
are  calculated  for  the  two-parameter  Weibull  distribution  using 
maximum  likelihood  estimating  equations.  Again,  if  there  are  cen- 
sored data,  an  iterative  procedure  is  used.  First,  Cohen's  vari- 
ance-covariance matrix  is  calculated.  Then,  from  Cohen's  vari- 
ance-covariance matrix,  Essenwanger's  variance-covariance  matrix 
is  calculated  and  the  parameters  for  Essenwanger's  two-parameter 
Weibull  distribution  are  calculated.  Using  the  PLOT  subroutine, 
the  Weibull  cumulative  and  probability  distribution  functions  are 
plotted.  If  there  is  more  than  one  sample  of  data.  Chi-square 
statistics  are  then  calculated  to  compare  each  sample  of  data  with 
each  previous  data  sample. 

4.  INPUT  PREPARATION.  The  deck  of  input  cards  follows  the  pro- 
gram EXECUTE  (@X0T)  card.  A FINISH  (OFIN)  card  is  the  last  card 
in  the  deck.  The  program  terminates  when  the  FINISH  card  is  read. 
There  are  two  types  of  input  deck:  one  for  data  sets  and  one  for 
use  when  data  is  to  be  generated  by  the  program.  Both  types  are 
described  below. 

a.  Oata  Set  Input.  Data  set  input  must  consist  of  four  card 
types.  The  data  are  grouped;  that  is,  if  two  or  more  data  values 
are  identical,  the  data  value  is  given  along  with  the  frequency, 
indicating  the  number  of  identical  data  values.  Each  type  of  in- 
put card  is  described  in  Table  1.  Card  types  3 and  4 are  repeated 
for  multiple  samples. 

b.  Data  Generation  Input.  When  the  data  are  to  be  randomly 
generated  by  the  program,  seven  card  types  must  be  used.  Each 
card  type  is  described  in  Table  2. 

5.  NUMERICAL  EXAMPLE,  a.  Problem  Description.  The  example 
problem  consists  of  68  completed  observations  and  32  censored  ob- 
servations. The  sample  was  progressively  censored  with  three  dif- 
ferent censor  times. 
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Table  1.  Description  of  Cards  for  Input  Data 


mi 

Format 

Program 

variable(s) 

Explanation 

1 

15 

I DATA 

IDATA=0,  read  data 

2 

15 

NR SAMP 

NRSAMP  = number  of  data  samples 

3 

215 

N.K 

N = number  of  uncensored  data 
groups 

K = number  of  censored  data 
groups 

4 

15. 

JFREQ(I). 

JFREQ!I)  = number  of  like  values 
in  the  i^^  group 

F10.5 

VALUE! I) 

VALUE! I)  = common  value  of  data 
in  the  i^"  group 

Table 

2.  Description  of  Cards  for  Randomly  Generated 

Data 

Card 

number 

Format 

Program 

variable(s) 

Explanation 

1 

15 

IDATA 

IDATA  = 1,  generate  data 

2a 

15 

IZZQ 

IZZQ  = number  of  random  numbers 
to  be  skipped 

3 

15 

I NO 

INO  = number  of  data  points  to 
be  generated 

4 

415 

NA,  NB, 

NC,  NS 

NA  = number  of  scale  parameters 

NB  = number  of  shape  parameters 

NC  = number  of  censor  percentages 
NS  = 0,  uniform  censoring 

= 1,  censor  upper  75  percent 

5 

6F10.2 

ALPHA! I) 

ALPHA! I)  = the  i^^  scale  parameter 

6 

6F10.2 

BETA! I) 

BETA! I)  = the  i^^  shape  parameter 

7 

6F10.2 

CEN!I) 

CEN!I)  = the  i^^  censor  percentage 

^This  card  eliminates  the  need  for  a new  random  number  seed  in 
the  uniform  random  number  generating  subroutine  by  skipping  IZZQ 
random  numbers. 
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b.  Program  Input.  The  program  input  consists  of  the  four  card 
types  listed  in  Table  1 and  four  systems  cards  required  by  the 
UNIVAC  1108  computer  (0RUN,  C^ASG,  OXgT  and  OFIM).  The  input  cards 
for  the  numerical  example  follow. 


Example 

aRUNf/TP  A106A#F1830A8397B*UNCLASSIFIC0t2  #2  00 
aAS6«A  06WEIBULL. 
aXOT  neWEIBULL.RUN 

c 

1 

47  3 

E.  177. 

1 246. 

1 252. 

1 269. 

1 283. 

1 294. 

5 331. 

1 367. 

1 379. 

1 386. 

1 411. 

1 423. 

1 441. 

1 488. 

1 502. 

1 508. 

1 519. 

1 531. 

1 542. 

1 550. 

1 553. 

1 568. 

1 583. 

1 589. 

1 601. 

1 613. 

1 621. 

5 682. 

1 772. 

1 806. 

1 820. 

1 340. 

1 854. 

1 872. 
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1 998. 

5 969. 
1 1033. 
1 1066. 
1 1088. 
1 1107. 

6 1184. 
1 1273. 
1 1309. 
1 1351. 
1 1374. 
1 1461. 
1 1494. 

10  246. 

15  742. 

7 1494. 

SPIN 


c.  Program  Output.  The  program  output  consists  of  the  uncen- 
sored  and  censored  data,  moment  and  maximum  likelihood  parameter 
estimates,  variance-covariance  matrix,  the  pdf,  the  cdf  and  other 
information.  Output  from  the  above  numerical  example  follows. 
Using  the  given  inputs,  the  shape  parameter  (S)  is  calculated  to 
be  1.91  and  the  scale  parameter  (ci),  993.24.  A shape  parameter  of 
1.91  gives  a positively  skewed  distribution  lying  between  the  ex- 
ponential distribution  (6  = 1)  and  the  normal  distribution  (6  = 
3.5).  On  the  cdf  graph,  the  ones  (1)  are  the  computed  cdf,  with 
the  twos  (2)  representing  the  observed  cdf  using  only  the  uncen- 
sored values,  i.e.,  ( E(uncensored  val ue)/(total  uncensored  + total 
censored  values).  Tne  threes  (3)  represent  the  ooserved  cdf  using 
all  values,  i.e.,  [ E(uncensored  values  + censored  values)]/(total 
uncensored  + total  censored). 


Example 


DATA  PITTING  OPTION  HAD  CHCSFN 

op.  GAMPLESr  1 
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